Load libraries.
library(Seurat)
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
Attaching SeuratObject
library(URD)
Loading required package: ggplot2
Loading required package: Matrix
Registered S3 method overwritten by 'gplots':
method from
reorder.factor gdata
library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.2 ──✔ tibble 3.1.8 ✔ dplyr 1.0.10
✔ tidyr 1.2.1 ✔ stringr 1.5.0
✔ readr 2.1.3 ✔ forcats 0.5.2
✔ purrr 1.0.1 ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::expand() masks Matrix::expand()
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
✖ tidyr::pack() masks Matrix::pack()
✖ tidyr::unpack() masks Matrix::unpack()
Load merged eb, lb, ehf.
load("/Volumes/Bud_SSD/20230122-urd/merged-seurat.Robj")
Grab count data and meta data.
count.merged <- as.matrix(merged.seurat@assays$RNA@counts)
Warning: sparse->dense coercion: allocating vector of size 6.2 GiB
meta.merged <- as.data.frame(merged.seurat@meta.data)
rm(merged.seurat)
Create urd.
merged.urd <- createURD(count.data = count.merged, meta = meta.merged, min.cells = 3, min.counts = 0)
2023-01-23 07:20:42: Filtering cells by number of genes.
2023-01-23 07:20:59: Filtering genes by number of cells.
2023-01-23 07:21:23: Filtering genes by number of counts across entire data.
2023-01-23 07:21:46: Filtering genes by maximum observed expression.
2023-01-23 07:22:11: Creating URD object.
2023-01-23 07:22:15: Determining normalization factors.
Warning: sparse->dense coercion: allocating vector of size 6.0 GiB2023-01-23 07:22:28: Normalizing and log-transforming the data.
as(<dgeMatrix>, "dgCMatrix") is deprecated since Matrix 1.5-0; do as(., "CsparseMatrix") instead
2023-01-23 07:23:12: Finishing setup of the URD object.
Warning: sparse->dense coercion: allocating vector of size 6.0 GiB2023-01-23 07:23:26: All done.
Add stage.
merged.urd@group.ids$stage <- as.character(merged.urd@meta[rownames(merged.urd@group.ids),"stage"])
Add embryo.
merged.urd@group.ids$embryo <- as.character(merged.urd@meta[rownames(merged.urd@group.ids),"embryo"])
#embryo <- sort(unique(merged.urd@group.ids$embryo))
embryo <- c("EB1", "EB2", "AVJ4-56", "AYQ8-11", "AYR4-57", "BBH6-17", "AXL7-21", "EHF1", "EHF2")
Calculate variable genes. Use only Flox embryos.
var.by.embryo <- lapply(c(1, 2, 3, 4, 5, 6, 7, 8, 9), function(n) {
findVariableGenes(merged.urd, cells.fit=cellsInCluster(merged.urd, "embryo", embryo[n]), set.object.var.genes=F, diffCV.cutoff=0.3, mean.min=.005, mean.max=100, main.use=paste0("Embryo ", embryo[n]), do.plot=T)
})
Add variable genes.
var.genes <- sort(unique(unlist(var.by.embryo)))
merged.urd@var.genes <- var.genes
Calculate PCA.
merged.urd <- calcPCA(merged.urd, mp.factor = 2)
[1] "2023-01-23 07:32:26: Centering and scaling data."
[1] "2023-01-23 07:32:34: Removing genes with no variation."
[1] "2023-01-23 07:32:36: Calculating PCA."
[1] "2023-01-23 07:37:28: Estimating significant PCs."
[1] "Marchenko-Pastur eigenvalue null upper bound: 1.40724420496983"
[1] "31 PCs have eigenvalues larger than 2 times null upper bound."
[1] "Storing 62 PCs."
pcSDPlot(merged.urd)
set.seed(19)
merged.urd <- calcTsne(object = merged.urd)
plotDim(merged.urd, "stage", plot.title = "tSNE: Stage")
Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as of ggplot2 3.3.4.
plotDim(merged.urd, "cluster.whole", plot.title = "tSNE: Cluster")
Save.
save(merged.urd, file = "merged-tsne-urd.Robj")
Load.
load("merged-tsne-urd.Robj")
Calculate diffusion map.
merged.urd <- calcDM(merged.urd, knn = sqrt(dim(merged.urd@meta)[1]), sigma=NULL)
[1] "destiny determined an optimal global sigma of 41.63"
Warning: You have 1603 genes. Consider passing e.g. n_pcs = 50 to speed up computation.as(<dsCMatrix>, "dgTMatrix") is deprecated since Matrix 1.5-0; do as(as(., "generalMatrix"), "TsparseMatrix") instead
plotDimArray(merged.urd, reduction.use = "dm", dims.to.plot = 1:8, outer.title = "Diffusion Map (Sigma NULL, 214 NNs): Stage", label="stage", plot.title="", legend=F)
Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as of ggplot2 3.3.4.
plotDim(merged.urd, "stage", transitions.plot = 10000, plot.title="Developmental stage (with transitions)")
plotDim(merged.urd, "cluster.whole", transitions.plot = 10000, plot.title="Developmental stage (with transitions)")
Add cluster and stage.
merged.urd@meta <- merged.urd@meta %>% unite("cluster.stage", cluster.whole:stage, remove = FALSE)
merged.urd@group.ids$cluster.stage <- as.character(merged.urd@meta[rownames(merged.urd@group.ids),"cluster.stage"])
Calculate pseudotime.
pseudotimePlotStabilityOverall(merged.urd)
plotDim(merged.urd, "pseudotime")
plotDists(merged.urd, "pseudotime", "stage", plot.title="Pseudotime by stage")
plotDists(merged.urd, "pseudotime", "embryo", plot.title="Pseudotime by embryo")
plotDists(merged.urd, "pseudotime", "genotype", plot.title="Pseudotime by genotype")
Save.
save(merged.urd, file = "merged-pseudotime-urd.Robj")
Extract pseudotime.
pseudo.data <- data.frame(rownames(merged.urd@pseudotime), merged.urd@meta$stage, merged.urd@meta$embryo, merged.urd@meta$genotype, merged.urd@meta$cluster.whole, merged.urd@meta$germ.whole, merged.urd@pseudotime$pseudotime)
save(pseudo.data, file = "pseudotime-urd.Robj")